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Abstract 

We simulate the Ising model on dynamical quadrangulations using a generalization of the flip 
move for triangulations with two aims: firstly, as a confirmation of universality for the KPZ/DDK 
exponents of the Ising phase transition, worthwhile in view of some recent surprises with other sorts 
of dynamical lattices; secondly, to investigate the transition of the Ising an£i-ferromagnet on a dy- 
namical loosely packed (bipartite) lattice. In the latter case we show that it is still possible to define 
a staggered magnetization and observe the anti-ferromagnetic analogue of the transition. 

Submitted to Physics Letters B. 



1 Theoretical Background 



The behaviour of various spin models on dynamical triangulations and their dual planar <j) 3 graphs 
has been extensively investigated both analytically |^| and numerically ||, || recently because of their 
interest for string theory and 2d gravity. The sum over triangulations or graphs in such models serves as a 
discretization of the sum over metrics in the continuum theories. In simulations the space of triangulations 
or their duals is generally sampled by performing local "flip" moves Js) , though more recently cluster-like 
moves for the geometry ("baby-universe surgery") have been proposed ||. Considering multiple copies 
of spin models, or alternatively large Q Potts models, has led to considerable progress in understanding 
the c > 1 regime in 2c? quantum gravity J?], ^| . 

The microcanonical partition function for an Ising model on some ensemble of graphs with n vertices 
is given by 

z n (p, # ) = E E ^ E + H E A c 1 ) 

G" { CT } \<ij> i J 

where G™- is the connectivity matrix for the graph G n . For the particular case of the Ising model on 
planar <f> 3 and </> 4 graphs the solution proceeds by noting the equivalence of the grand canonical partition 
function arising from equ. (Q) and the free energy of an N x N two matrix model . In the (f> 4 case in 
zero external field this is 

F 4 (c, g,H) = — log ( J d N Mid M 2 exp -tr- ( -± + -± - cM x M 2 + -M* + -Af 2 4 J J (2) 
and for <j) 3 graphs 

(c, g, H) = ± log (| d N2 Ml d N2 M 2 exp -trj + ^ - cM x M 2 + hi 3 + ±M 3 ^j ) (3) 

where c = exp(— 2/3) introduces the temperature of the spins into the theory. In the planar limit N — > oo 
it is possible to solve the model exactly using orthogonal polynomial or saddle point methods and derive 
the full set of critical exponents, which are in agreement with those calculated from continuum conformal 
field theory §. 

It is expected that the Ising model on any dynamical polygonization or its dual will display the 
KPZ/DDK critical behaviour so long as further couplings are not tuned to give a multicritical point. 
This has been demonstrated both analytically and numerically for various variants of 4> 3 graphs where 
self-energy diagrams and tadpoles are either excluded or not. Although the critical temperature changes, 
the critical exponents are unaffected. There have, however, been some recent surprises with other dis- 
cretizations: with 2d Regge calculus the Onsager exponents were observed ||; and with apparently "flat" 
models the KPZ/DDK exponents were measured ]Toj . This increases the incentive for performing an 
explicit check on the behaviour of various possible discretizations. 

It is possible to write down explicitly a matrix model for the Ising model on a dynamical triangulation 
(ie with the spins residing at the vertices of the triangle, rather than in the centre as is effectively the 
case for tfi 3 graphs). The action of this model is given by 

U DTR s = -tr(- -±——S 2 + - — L -D 2 + S 3 /3 + SD 2 ) (4) 

g \2cosh(p)(l + c*J 2 cosh(p)(l — c*J / 

where c* = (1 — c)/(l + c) is the transformation induced in c by the standard duality transformation on (3 
and S and D are still N x N hermitian matrices representing edges with the same and different spins at 
the end respectively. This model may be transformed into a form equivalent to the "O(l)" representation 
of the Ising model on <fi 3 graphs which is discussed in 

^M^+C-^+ ^-y^'M (5) 
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and shown to have Ising critical behaviour at the dual of the <p 3 model transition temperature. Working 
directly on the triangulations with Udtrs thus provides an explicit confirmation of duality and a further 
check on the universality of the exponents. 
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We can carry out a strong coupling expansion whatever the lattice our Ising model inhabits 

Z v ((3, 0) = 2 y (cosh(/3)) B ]T (tanh(/3))« (6) 

G v loops 

where the loop sum is restricted to edges that are only traversed once. This makes it clear that we should 
only expect to see an anti-ferromagnetic transition in the Ising model on lattices where only even loops are 
possible, because only here do we have a (3 — > —(3 symmetry. Thus, for instance, the honeycomb lattice 
in 2d displays an antiferromagnetic transition whereas the triangular lattice does not. It is not entirely 
trivial to conclude that such transitions would survive a coupling to 2c? gravity, which in effect introduces 
annealed connectivity disorder into the lattices. Both the analogue of the honeycomb lattice - dynamical 
cj> 3 graphs, and the triangular lattice - dynamical triangulations, fail to display an antiferromagnetic 
transition because both odd and even loops are present. 

Graphs with only even loops can be constructed in general by considering complex matrices ]l3[ , but 
we take a slightly different tack here and consider the Ising model on a dynamical quadrangulation, where 
the action is given by 

U DQRS = jtr Qs 2 + \d 2 + |s 4 + |i? 4 + 9 -l{SDSD + 2S 2 D 2 ^ (7) 

with 

51 = cosh 2 (/3)(l + c*) 2 

52 = cosh 2 (/3)(l-c*) 2 

53 = cosh 2 (/3)(l + c*)(l-c*). (8) 

This model can be shown to have a transition at c = 3/5 by transforming the action back to a form 
resembling that for the Ising model on <ft 4 graphs. With some scalings we can rewrite equ.(^) as 

U DQRS = £tr Q* 2 + \ D 2 + —±^S* + + ^^— iSDSD + 2S 2 D 2 )) . 

(9) 

The model of equ.(g) is recast by setting A = (Mi + M 2 )/V2, B = (Mi - M 2 )/\/2 and then rescaling 
A -> V2A/y/l~~i, B -> v^S/yTTc, g -> 2g. This gives 

[/, 4 = —tr ( -A 2 + -B 2 + —J— A 4 + ; 1 R 4 + L -(ABAB + 2A 2 B 2 )] . (10) 

5 V2 2 4(1 -c) 2 4(1 + c) 2 2(l-c)(l + c) V ') y ' 

which is identical in form to equ. ([)]). Given that £7^4 produces an Ising transition at c = 1/4, comparison 
with equ.(^|) shows that Udqrs produces an Ising transition at c* = 1/4 (ie c = 3/5), the quoted result. 
Furthermore, the possibility of an antiferromagnetic transition is manifested in equ.(j^) too as (3 — > —j3 
(ie c* — > — c*) is obviously a symmetry of the action upon exchange of S and D. This is not the case for 
the original 4 action, even in its modified form in equ.(|Tc|), as c — > 1/c when (3 — ► — [3. 

The matrix model considerations would thus suggest that the Ising model on dynamical quadrangula- 
tions displays an antiferromagnetic transition with the same critical exponents as the ferromagnetic case. 
This is not a totally vacuous observation, as the antiferromagnetic transition essentially operates on two 
decoupled sets of Ising spins, so it is not entirely obvious that one should see the exponents for a one Ising 
model transition coupled to gravity rather than those for two Ising models. To check our observations 
with a simulation requires the definition of a flip move in the space of quadrangulations, that possesses 
all the nice ergodicity properties of the standard triangle flip moves. In | fll) we baldly presented one 
candidate move, shown again in Fig.l, without any discussion of its derivation or ergodicity properties. 
We now remedy this lacuna. 

The easiest approach to take is to break our quadrangulation into a triangulation and reconstruct 
the possible moves from that. We shall consider undecorated (no spin) quadrangulations first. Not every 
triangulation can be paired off to give a quadrangulation, we need one that admits a perfect matching. 
A matrix model to do just this (in the dual </> 4 /^ 3 picture) was discussed in JlJ] 

Z = J D" 2 J> exp- Q^ 2 - ^e 1 " 2 -^ 4 ) = J D N ' > D N \ exp- Q^ 2 + \u 2 -J= e -^ 2 U ) 

(H) 
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where both <f) and u are N x N hermitian matrices. Integrating out u on the right hand side gives back 
the single matrix model on the left hand side. If we go back to the direct lattice picture of triangles and 
squares the effect of the integration in the above dual lattices is to pair off all triangles into squares. If we 
look at the dual <j) 4 diagram to two adjoining squares it can be broken down, or perhaps more correctly 
exploded, in four possible ways with the insertion of u propagators as shown in Figs. 2. The standard 
<j) 3 flip moves can be performed subject to the constraint that every vertex must retain one u and two </> 
edges, before collapsing the u propagators to obtain the final (j) 4 diagram. Performing the allowed flips 
on Figs. 2c, d leave the final (f> 4 diagram unchanged and the results of the allowed flips on Figs2a,b are 
shown in Fig. 3a and Fig. 3b. Just as in the </> 3 flip move both labellings of the central segment appear 
with equal probability in both the diagrams, and just as for the (j) 3 flip the distinction disappears when 
one returns to the direct lattice. 

When decorating the lattice with spins it is best to think of the spin at the exploded vertices living 
on the u propagator and getting carried with it in any flip moves. This avoids any potential conflict that 
might arise when collapsing the diagram after performing the flips if we actually created new spins on the 
split vertices. The net result of all of this is that the flip moves deduced by the procedure of exploding the 
<j> 4 diagram into (j) 3 diagrams and carrying out (j) 3 flips before reconstructing the <j> 4 diagram are precisely 
those given in Fig.l. Both the moves appear with equal probability. As we can trace the moves back to 
the <j) 3 picture (albeit with two sorts of propagators) the ergodicity is assured. 

A final subtlety to consider before launching into the simulations proper is the choice of class of 
quadrangulation. We restrict ourselves for reasons of computational simplicity to non-degenerate quad- 
rangulations in which we do not allow such nasties as squares touching on more than one side. In the 
dual <j) 4 picture this corresponds to excluding bubbles on propagators, the "setting-sun" diagram and the 
one loop contribution to the <f) 4 vertex. Recent simulations on (j) 3 graphs have shown, rather counterintu- 
itively, that including the (j) 3 counterparts of such diagrams actually improves the scaling behaviour and 
we hope to return to this issue for quadrangulations in further work. 

2 Simulations 

We simulate graphs of size N = 100, 225, 400, 900 and 2500 vertices. In all cases the starting configuration 
was taken to be a regularly quadrangulated torus for simplicity (e.g. 2500 = 50 x 50). The topology 
of the manifold on which the spin model lives is not expected to change the KPZ/DDK exponents for 
the spin model transition, although j string, which we do not measure here, will of course depend on the 
topology. We conducted 10,000 equilibration sweeps and 50,000 measurement sweeps for most of the 
lattice sizes at each (3 value simulated, with more sweeps on the larger graphs and closer to the phase 
transition. For the ferromagnetic model we updated the spins with the Wolff cluster algorithm, ensuring 
that on average every spin was updated per sweep. The number of flip updates per sweep was chosen 
to be equal to N, a rule of thumb that works well for <j) 3 graphs and triangulations, and seems to do so 
here too. Varying the number of flip updates by a factor of two in either direction had little effect on the 
measured quantities. Although it is possible to construct a cluster algorithm for the antiferromagnetic 
model as well, we contented ourselves with an old-fashioned Metropolis update in this case, again for 
simplicity as very high accuracy was not a primary requirement. 

We measured all the standard thermodynamic quantities for the model: the energy E, specific heat C, 
magnetization M, susceptibility \ and various cumulants. A short word is in order on the magnetization 
and energy measurements for the antiferromagnetic model. The normalization of the energy (ie the 
constant term in the hamiltonian) is chosen so that the energy runs between 2 and 4 for both the 
ferromagnct and antifcrromagnct. In Fig. 4 we plot the energy in both models. It is clearly essentially 
identical, especially when the less efficient algorithm employed for the antiferromagnet is taken into 
account. This provides the first confirmation of the matrix model result that the antiferromagnetic 
transition should be identical to the ferromagnetic transition on dynamical quadrangulations. We have 
to work slightly harder to demonstrate the equality of the magnetizations as the standard magnetization 
will clearly stay zero in any antiferromagnetic transition. It is, however, still possible to easily define a 
staggered magnetization on the dynamical quadrangulation because of our choice of starting configuration. 
We can assign an initial sign or parity to each vertex in a "checkerboard" pattern consistently on the 
initial regularly quadrangulated torus. The flip move of Fig.l maintains the parity of flipped links as is 
clear from looking at the solid and open dots in the diagram, so we can define a staggered magnetization 
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as 

M =^E(- 1 ) i+ ^ (12) 

where i,j labels the original horizontal and vertical coordinates of the spin in the starting torus. Al- 
though our choice of a regular torus as starting point makes the definition of the staggered magnetization 
particularly easy, it should be stressed that any loosely packed lattice could be checkerboarded in such a 
manner and that the corresponding flip moves would not disrupt these assignments or the definition of 
the staggered magnetization. In Fig. 5 we plot the standard magnetization for the ferromagnet and the 
staggered magnetization for the antiferromagnet, which are again very nearly equal until the inefficiency 
of the metropolis algorithm begins to bite in the low temperature phase. Note the large error bars in this 
case; all the other error bars are smaller than the symbols. 

Having convinced ourselves of the qualitative similarity of the ferromagnetic and antiferromagnetic 
transitions from these figures we now proceed to a more quantitative scaling analysis. We commence by 
looking at the Binder's cumulant for the magnetization 

whose crossing point for different lattice sizes determines (3 C . The data for both the ferromagnet and 
antiferromagnet give quite clean crossings that are consistent for the various different size pairings, in 
contrast to our </> 3 simulations in Q] where we had to resort to a slightly more devious approach to 
extract a reliable estimate for (3 C . We find C — 0.398(2) for the ferromagnet and j3 c — —0.40(2) for 
the antiferromagnet, where the larger error in the latter is due to the poorer statistics of the metropolis 
algorithm employed. Note that we have not found j3 c — — l/21og(3/5) ~ 0.255 ... as given by the matrix 
model because the matrix model calculation implicitly assumes the inclusion of the degenerate squares 
that we have expunged. It would be an interesting exercise to calculate the expected (3 C for our variant 
analytically to compare with the measurements above. 

We can extract the combination vD from Um in various ways. One possibility is to look at the 
maximum slope of Um , which is expected to scale as 



max 



fdU A 



V d0 



N x ' vD . (14) 



For the ferromagnet we find that 1/vD = 0.31(5) or vD = 3.2(5), where the rather large error bar is 
accounted for by the effect of the smaller lattices. Fitting to only the larger lattice sizes gives both a 
slightly smaller estimate and a smaller error. If we accept the likelyhood of finding KPZ/DDK rather 
than Onsager exponents the hyperscaling relation a = 2 — vD would give us vD = 3 rather than vD = 2. 
The above value is clearly in better agreement with the KPZ/DDK result than the Onsager result. An 
alternative is to look at the scaling of the slope at the crossover point, which gives slightly a slightly 
better fit of 1/vD = 0.31(3) or vD = 3.2(3), consistent with the above values. The results for the 
antiferromagnet are very similar, again with larger errors due to the poorer statistics. Given vD and P c 
we can now proceed to a finite size scaling analysis of some of the other exponents. In the discussion 
which follows, the results are for the ferromagnet unless explicitly indicated otherwise. 

Although a fit to the finite size scaling relation for the specific heat is not particularly illuminating 
because of the extra adjustable constant A that is present 

C ~ A + BN a ' vD (15) 

it is consistent with the values of vD coming from the cumulant. The specific heat for the ferromagnet 
plotted in Fig. 6 (that for the antiferromagnet is identical) clearly shows the expected cusp, the weak 
growth of the peak for small lattice sizes dying away as the size is increased. Working in lexicographic 
order, we now turn our attention to the exponent 

M ~ M N- p/vD (16) 

which appears to give rather poor fits on modestly sized 4> 3 lattices [0 , although the use of a triangulation 
H or the inclusion of tadpoles and self-energies greatly improves matters JtJ . With quadrangulations we 
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find that we have maintained our traditionally atrocious fits to 0, finding fi/vD = 0.075(2), which 
translates to (3 = 0.24(3) rather than the KPZ/DDK value of 0.5. The fits to 7 from the susceptibility x 

X * X0 N^ vD (17) 

are much more satisfactory, however. We find j/vD = 0.72(2), which using our best fitted value of 
vD = 3.2(3) gives the estimate of 7 = 2.3(3) in quite good agreement with the KPZ/DDK value of 2. 
We can see clearly the sharp divergence in x m Fig. 7. A similar story holds for the antiferromagnetic 
simulations: a good fit to the KPZ/DDK value for a/vD is still possible because of the constant term, 
we find fi/vD = 0.078(2) - again too small, and 7/1/Z) = 0.73(7). 

The lattices we simulate are rather modestly sized, but it is still possible to do a direct fit to the 
dependence of C, M and x on the reduced temperature t = \(3 — c \/Pc on the largest lattice (N = 2500) 
as a consistency check. If we fit 

C~A' + B't~ a (18) 
we again get a good value a = —1.1(1) because of the extra constant A' . The fit to j3 

M ~ M'tf 13 (19) 

is even more mediocre than the finite size scaling fit for f3/vD, giving a wide range of values with similar 
goodness of fit depending on the points deleted from the fit. The direct fit to 7 

X^Xot~\ (20) 

on the other hand, is very good: 7 = 2.07(2). 

For completeness we list our fitted exponents and critical temperature for the ferromagnetic simula- 
tions below along with the analytical values for the critical exponents 





a 


P 


7 


vD 


ctjvd 


0/ud 


j/ud 




measured 


-1.1(1) 




2.07(2) 


3.2(3) 


-0.35(2) 


0.075(2) 


0.72(2) 


0.398(2) 


theory 


-1 


1/2 


2 


3 


-1/3 


1/6 


2/3 





Table 1: Measured and Theoretical values of a,/3, 7, vD. 
(We have not included the very poor direct fit to (3.) 

We close with a few observations on the efficiency of the Monte-Carlo algorithms. Previous work on 
triangulations and 4> 3 graphs has shown that although cluster algorithms greatly reduce critical slowing 
down for energy observables the magnetization still has quite large values of the dynamical critical 
exponent z M , defined through the scaling of the appropriate autocorrelation time r 

t~N z ' d . (21) 

For the ferromagnetic simulations with the Wolff cluster algorithm we find zm/D = 0.39(3) for the 
magnetization, and ze/D — 0.23(4) for the energy. These values are unlikely to be particularly accurate 
because of the small lattice sizes, but the value for the energy is unusually large compared with previous 4> 3 
and triangulation simulations. A glance at the metropolis simulations for the antiferromagnet shows that 
the cluster algorithm autocorrelation times are an order of magnitude better for the magnetization and 
a factor of two or three better for the energy on the N = 2500 graphs, so they still offer an improvement 
over a local algorithm in spite of the anomalously large Ze- 

It was suggested in || that the bottlenecks that appear in the baby universe structure of the 4> 3 graphs 
and triangulations impeded cluster formation and was responsible for this relatively poor performance. 
Although we have not investigated the baby-universe distribution of our quadrangulations there is no 
reason to suppose it is significantly different from triangulations as we are clearly simulating the same 
continuum theory if the Ising results are to be believed. We would therefore expect the same bottlenecks 
to be present on our quadrangulations, causing the same obstacles to cluster growth. By the same token, 
it would be equally easy to perform baby universe surgery Q on a quadrangulation to help alleviate the 
problem, the only difference from a triangulation being that the minimum neck is now of length four 
rather than three. 
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Other geometrical features such as the flip acceptance as a function of temperature and the cluster 
distributions look qualitatively similar to earlier tfi 3 and triangulation simulations. As an additional 
measure to improve our measurements of the magnetization critical exponent it would be a relatively easy 
matter to include degenerate quadrangulations in the class of graphs we simulate in order to improve the 
convergence to the continuum limit. 

3 Conclusions 

Modulo the usual rather poor fits for the magnetization exponent (3 the measured critical exponents for 
"square gravity" are consistent with the KPZ/DDK values, thus providing further explicit numerical 
support for the notion of universality on dynamical lattices as far as the ferromagnetic transition is 
concerned. As a by-product of the simulations we have also found the critical temperature for the non- 
degenerate quadrangulations that we considered here (3 C — 0.398(2). This translates to the dual value 
(3 C ~ 0.486 on non-degenerate (j) 4 graphs, to be compared with f3 c ~ 0.693 for (j) 4 graphs in which the 
various bubbles etc. are allowed. 

It is a relatively simple matter to simulate the dual <fi 4 model as well, but here of course we would 
no longer find an antifcrromagnetic transition. We need to consider loosely packed lattices on which 
Neel order is possible to observe the antifcrromagnetic transition on fixed lattices, which remains the 
case for dynamical lattices. On dynamical lattices we are summing over some class of loosely packed 
lattices, which we implement in the simulations with a compatible flip move. From the matrix model 
and the simulations it is clear that the antifcrromagnetic transition survives its coupling to 2d gravity 
on dynamical quadrangulations and retains the same exponents as the ferromagnetic transition. The 
antifcrromagnet thus provides an amusing example of Neel order surviving on what is essentially a fluid 
random lattice. 

Extending the approach here to investigate other polygonizations is also not difficult. Exploding 
higher order vertices into (j) 3 vertices in such models allows a derivation of the appropriate flip moves and 
weighting factors. 
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Figure Captions 

Fig.l The two possible flip moves on adjacent squares in DQRS both preserve the "checkerboard" order 
of the ground state. Up spins are shown as solid dots and down spins as open dots. 

Figs. 2 The possible ways of "exploding" the (f> 4 diagram dual to two adjacent squares. The u propagators 
are represented as dotted lines and the (j> propagators as solid lines. 

Figs. 3 The resulting <j> flip moves once the diagram has been pieced together again after carrying out 
all the possible flips in the exploded model. Both orientations of the central link occur with equal 
probability in both diagrams. Going back to the direct lattice gives the move depicted in Fig.l 

Fig. 4 The energy for both the ferromagnetic and antiferromagnetic models on the largest lattices simu- 
lated [N = 2500). 

Fig. 5 The magnetization for both the ferromagnetic and antiferromagnetic models on the largest lattices. 
Fig. 6 The specific heat C for the ferromagnet. 
Fig. 7 The susceptibility \ for the ferromagnet. 
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